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Abstract. In this paper we report on the results of a numerical study of the nonlinear time-dependent 
Bardeen-Cooper-Schrieffer (BCS) equations, often also denoted as Bogoliubov-de-Gennes (BdG) equa¬ 
tions, for a one-dimensional system of fermions with contact interaction. We show that, even above the 
critical temperature, the full equations and their linear approximation give rise to completely different 
evolutions. In contrast to its linearization, the full nonlinear equation does not show any diffusive behavior 
in the order parameter. This means that the order parameter does not follow a Ginzburg-Landau-type 
of equation, in accordance with a recent theoretical result in [T]. We include a full description on the 
numerical implementation of the partial differential BGS/BdG equations. 
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1 Introduction 

When Bardeen, Cooper and Schrieffer, shortly BCS, pub¬ 
lished one of the most famous papers in physics in 1957 [5], 
giving the first microscopic explanation for superconduc¬ 
tivity, a phenomenological theory for the phenomenon had 
already been around. Systems close to the critical tem¬ 
perature were described with the help of a macroscopic 
phase-transition parameter introduced by Ginzburg and 
Landau in 1950 [5]. Their theory was the first one to al¬ 
low for the description of the spatial dependence of the 
superconductivity inside superconducting alloys and the 
first with which to explain type-II superconductors and 
the hexogonally shaped penetrations by magnetic flux. 

As the Ginzburg-Landau theory yields reliable results 
on the large scale, soon the question arose as to whether 
this model can be understood as a macroscopic limit of 
BGS theory for systems close to the critical temperature. 
Gorkov gave a positive answer to this question for the sta¬ 
tionary case shortly after the publication of BGS, see [6]. 
A rigorous mathematical proof of the convergence was 
achieved some years ago [7]. 

But what remains unclear and controversial up to this 
day, in particular in terms of a rigorous derivation, is the 
question whether the time evolution of superconducting 
systems close to the critical temperature are governed 
by a Ginzburg-Landau type of equation. After first at¬ 
tempts for a derivation of the macroscopic limit had been 
presented Gorkov and Eliasberg pointed out 

that a nonlinear equation could only be valid in a gapless 
regime m- Still, in the authors made a case 

for a time-dependent Ginzburg-Landau equation for su¬ 
perfluid gases at temperatures slightly above the critical 


one. The argument is based on the assumption that the 
nonlinear terms in the BGS/BdG equations only lead to 
small perturbations but do not quantitatively change the 
system’s behavior. In more detail, this would mean that 
the projection of the Cooper pair density onto the cen¬ 
ter of mass direction is governed by a nonlinear dispersive 
equation. However, it has been argued recently in [T] that 
for a translation invariant homogeneous system close to 
equilibrium, the full BGS/BdG equations and their lin¬ 
earization do not yield the same behavior at temperatures 
close to the critical one. In particular, dissipative behavior 
can only be expected for the linear approximation of the 
BCS equations but not for the full equations. 

With our work, we demonstrate this result by means of 
a thorough numerical study of the long-term evolution of 
the BCS equations and their linearization for spatially ho¬ 
mogeneous systems close to equilibrium at temperatures 
slightly above the critical one. For decreasing values of 
the parameter h, defined via T = (1 -|- h?)Tc^ we evolve 
the full and the linear system over a long time span and 
track the behavior of the Cooper pair density and the or¬ 
der parameter. For each values of the small parameter h, 
we find clear differences between the full equation and its 
linearization. Additionally, we see that the full BCS /BdG 
equations yield oscillations in the order parameter about 
a constant value. Such a behavior has long been predicted 
for and already been observed in out-of-equilibrium sys¬ 
tems, see, e.g., [3SJIM1I37]. Although the focus of our study 
is not on oscillations in particular but rather on the long¬ 
term behavior of the equations in general, it is interesting 
that we can replicate such oscillations for systems close to 
equilibrium. 
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In the realm of numerical analysis, the treatment of 
quantum dynamical systems has been of huge interest for 
many decades (see [21] for an extensive overview). Vari¬ 
ous evolution schemes for the linear Schrddinger equation 
in varying settings have been proposed, see, e.g.. 

Nonlinear Schrddinger equations such as the 
Gross-Pitaevskii equation and equations arising from the 
Hartree and Hartree-Fock approximation of the quantum 
state have also been devoted attention to, see, e.g., Eliiii 
[3011^ and [^1^ . Regarding the BCS regime, the sta¬ 
tionary equations have been treated numerically in [33] 
and the time-dependent BCS/BdG equations have been 
considered from an analytical perspective in [34]. But, 
the above-mentioned studies of the out-of-equilibrium dy¬ 
namics of the BCS equations notwithstand¬ 

ing, to the best of our knowledge the coupled nonlinear 
time-dependent BCS equations have not been paid much 
attention to from a numerical point of view. Therefore, 
we come up with a reliable integration algorithm for the 
evolution of the system. 

The paper is organized as follows: First, we introduce 
the system we are considering and the physical background 
in Section |2| This is followed by a brief summary of the 
theoretical results of [T] in Section |3| Then, we present 
our numerical results for the linear and the full equation 
in Section |4| Finally, we summarize our main results in 
Section |5| A detailed discussion of the initial setup of the 
system and the numerical implementation is provided in 
the Appendix Sections [A] and jB] respectively. 


where the entropy S is defined as 

S{r) :=- [ Trc 2 (F(p) logF(p)) dp. (3) 
Jr 

The evolutions of a and 7 are given by the time-dependent 
BCS equations which are also known as Bogoliubov-de- 
Gennes equations |39j . In momentum space they can be 
written conveniently in the self-consistent form 

iFtip) = [HrAp),rt{p)]. (4) 

Here, the subscript t indicates the time-dependence and 
the Hamiltonian Hpt (p) is defined as 


HrAp) 


( ZJ^ 2[V*dt](p)\ 
\2[V*at](j)) p-p^ / 


(5) 


with * denoting the convolution. Calculating the upper- 
left and upper-right entries of the matrix-valued equa¬ 
tion ®, we arrive at the system of coupled nonlinear 
equations 


^lt{p) 


2 


{V * at){p)at{p) - {V * at){p)&t{p) , 


( 6 ) 


latip) = 2(p^ - pi)at(p) + 2{V * at){p) - 4(V + at){p)^t{p)- 

(7) 


2.2 Contact interactions 


2 Physical and mathematical background 

2.1 Energy functional and BCS equations 

In mathematical terms, BCS theory is a special case of 
a generalized Hartree-Fock variational principle, itself de¬ 
scribed by Bogoliubov-theory, for the density operators 
7 : H "H and a : ^ "H acting on the considered 

Hilbert space T-L. Those matrices are put together to form 
the two-by-two operator-valued matrix 

see, e.g., [38] for an introduction. The entries of the matrix 
can be represented by means of their momentum distribu¬ 
tion 7 (fc) = (a^Ofe) and the pair density a{k) = (akCL-k), 
determining the Cooper pair wave-function via Fourier 
transform as a{x — y) = (27r)~^^^ f We 

suppress spin in our notation; the pair density d is as¬ 
sumed, for simplicity, to be a spin singlet. For a one¬ 
dimensional translation invariant system of fermions at 
temperature T interacting via a potential V, the BCS 
pressure functional per unit volume is given by 

^T{r) = [ {p'^ - p)l{p)<^P + [ |Q:(a;)pf^(a;)dx - TS'(F), 

Jr Jr 

( 2 ) 


In this paper we concentrate on attractive contact inter¬ 
actions, i.e., potentials of the form 

V{x) = —a5{x), a > 0, ( 8 ) 

which lead to exactly solvable systems in the stationary 
case. Not only is such a potential the most interesting one 
from a physical model point-of-view but also does it allow 
us to implement the terms including a convolution in the 
equations of motion conveniently as we will illustrate in 
the numerics Section IS 

2.3 Initial values 

In this work we consider initial data which, in the sta¬ 
tionary case, could be described by the Ginzburg-Landau 
energy functional for temperatures T close to the critical 
temperature T^, i.e., T = Tc + h? for a small parameter 
h € R. For temperatures above Tc, the free energy is mini¬ 
mized by the so-called normal state for which = 0 , 
7 ^ = i/i-i-exp((p^-M)/T). For initial data /q to be within the 
range of Ginzburg-Landau, they have to satisfy 

TT{ro)-W^)<o{hA. (9) 

This condition can be complied with by choosing 

° 1 + e^^o/T 


( 10 ) 
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with 


Hao 


'p^ ^ -^0 \ 
-Ao 


( 11 ) 


where Aq is a small parameter of the order of h, see, e.g., 
[Tj- Calculating the right-hand side of the matrix equa¬ 
tion gives 


Co = 


7o ap \ 
ao 1-70/ 


where % and dp take the special form 


( 12 ) 




tanh 

~7W^W+W 


do 


tanh 


f ^( p 2_^)2 + | 4|2 

I 2T 


./w^w+w 


(13) 

(14) 


In our simulations we choose a temperature which is slightly 
above the critical temperature for the setting under con¬ 
sideration and set the initial value for the gap parameter 
Aq to a non-vanishing value. We explain how to obtain 
the critical temperature for our setting and how to find 
physically reasonable initial values for A in the Appendix 
Section s 


2.4 Ginzburg—Landau and macroscopic parameter 


For the stationary case it is well known that the Ginzburg- 
Landau theory emerges as the macroscopic limit of the 
BCS theory. To be more specific, define |a*) as the trans¬ 
lation invariant minimizer of the BCS functional which, in 
case of the contact interaction ([5]) , can be calculated via 


a*{p) 




2T 




(15) 


Then, for the Cooper pair density |a) corresponding to the 
non-translation invariant minimizer of Tt, the quantity 


:= 


1 

h 


(ala) 


(16) 


with some appropriate parameters cgl,i and cgl,2, see, 
e.g., [18] and [20l Eq. (18)]. The parameter cgl.i depends 
on i'r-Tc)/(h'^Tc). Crucially, cgl,i has the same sign as 
{T — Tc). Thus, the TDGL equation is dissipative for tem¬ 
peratures above Tc by definition. This implies that if ipt 
could be described by the TDGL for small h it should de¬ 
cay over time. However, we will demonstrate in Section |4| 
that this is not the case, at least for the full non-linear 
equation. The same conclusion has been reached by an 
analytical investigation recently as we will outline in Sec¬ 
tion O 


2.5 The linear approximation 

Let us decompose the particle density as 

7t=7o+??t- (19) 

For states satisfying ([Hj) pt appears to depend quadrati- 
cally on at, see, e.g., [U Eq. 11], and it seems legitimate 
to approximate the full equation by its linearization 

iat{p) = 2{p^ - p)at{p) + 2{V * dt)(p)(l - 2%{p)). 

( 20 ) 

However, close to the Fermi-surface the quantity pt is 
not small but the dominant part in the non-linear evo¬ 
lution. Consequently, the full BCS equations (Hj)-© and 
the linearization (EOj) give rise to very different evolu¬ 
tions. Namely, Eq. (l20ll yields a dissipative behavior in 
'ipt whereas the full equations do not as is shown formally 
in [T] and as we confirm by our numerical experiments be¬ 
low. Let us briefly summarize the results of [I] in the next 
Section. 


3 Recent mathematical results 

The BCS time-evolution dH) is studied analytically in [T]. 
Based on the work |7] the authors prove in [T] Theorem 
1] that I'l/'tj does not vanish for any times. More precisely, 
it is shown in a very general setting that, if the initial 
state Fq is close to the energy of the normal state, i.e., 
^riFo) — Frir^) < 0{h'^), then the corresponding tpt 
satisfies 


is an approximate solution of the stationary Ginzburg- 
Landau equation, see, e.g., m- This told, if there were 
an analogous relation between the time-dependent BCS 
and the CL equations, the order parameter 

ipt ■■= ^{ a *\ at ) (17) 

should, close to Tc, approximately satisfy a conventional 
time-dependent Ginzburg-Landau (TDGL) equation. In 
the spatially homogeneous case we are studying in this 
work, the conventional TDGL equation takes the form 

(18) 


( 21 ) 

for an appropriate constant C independent of h. 

On the other hand, it is shown in [1] that the solution of 
the linearized equation (1^ tends to 0 exponentially fast 
compared to the system’s time scale of In detail, using 
strategies from perturbation theory, it can be derived that 

l^tl (22) 

holds, where A is a resonance of order i//i^ which emerges 
from the zero-eigenvalue at T = Tc of the linear operator 

O = {k^ — p) tanh”^ ( ^2T^ ) 


Ipt = -CGL.lV’t - CGL,2|^/’tPV’ti 
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The combination of the bounds m and (EH) shows 
clearly that the non-dissipative behavior of '0t is a purely 
non-linear effect which takes place solely in a tiny neigh¬ 
borhood of the Fermi surface. 

Furthermore, using the methods of [T], it is straight¬ 
forward to derive the following bound on the derivative 

\^t\=0{yh). (23) 

In other words, although the solution |'0t | tends to the con¬ 
stant 1 ^ 0 1 in the limit h —0, its derivate might well oscil¬ 
late more and more -in line with according predictions for 
systems which are suddenly perturbed out of equilibrium 
l|35]b These findings are well reproduced in our numerical 
experiments as we show now. 

4 Simulations 



Fig. 1. The gap Z\o as a function of the semiclassical parameter 
h in semilogarithmic scale. 


In this work we are interested in a qualitative study of the 
differences between the full BCS/BdG equations and their 
linearization. Thus, without loss of generality, we can work 
in dimensionless units and set the constant a of the contact 
interaction and the chemical potential // to a = I and 
/r = I, respectively. The initial data for the simulations 
are obtained as outlined in the Appendix Section [A) For 
this, we approximate the integrals in Eqs. (OSl) and (1^51) by 
the sum over the discrete momenta we take into account. 
For the sake of reproducibility we add the thus-obtained 
values for and Aq to the results of our simulations. 
For more details on the discretization of the equations 
under consideration we refer the interested reader to the 
Appendix Section IbI 

4.1 Gap as a function of h 

In order not to have to calculate the initial value of the 
gap parameter which depends on the crucial parameter 
h at the start of each evoluation again, we calculate Aq 
with the procedure outlined in the Appendix Section 
for various h once. The interesting result is illustrated in 
Fig. □ where we can see that Aq depends more or less 
linearly on the crucial parameter. 

Finally, with both and Aq at hand, we are able to 
present the results of the simulations. Doing so, we take 
into account that at temperatures T = Tc + h?, physically 
interesting dynamics are expected to occur on a time-scale 
of Oy/hy. Therefore, we always set tend = ^end = 

2/^2 in the following. 


4.2 Results for h = 1/4 

We plot the scaled i^-norm of a, which in the discrete 
setting is given by the sum over the K discrete momenta 
as 




K/2-1 

E 




(24) 





The physical parameters are Tc = 0.19 and Ao — 0.29. 


as well as the modulus of the interesting macroscopic pa¬ 
rameter Tpt introduced in Eq. dnD- We plot the results 
for both the BCS equation © and its linear approxima¬ 
tion (l^ni) . see Fig. [5] and Fig. [3] For both quantities, the 
linear equation leads to exponential decay. The full equa¬ 
tion, in contrast, coincides with the linear approximation 
only for a short period after which both ||at|| and \tjjt\ 
grow again. 


4.3 Results for h = i/s 

Here, too, we consider the scaled norm of the Cooper pair 
density and the modulus of the parameter The results 
are shown in Fig 2] and Fig. [5] Again, the linear evolu¬ 
tion equation is clearly diffusive while the full equation 
yields a similar behavior only for very small times. After 
a short decline in the beginning of the simulation, ||at|| 
and liptl seem to oscillate. Similar oscillations have been 
predicted by [35] and observed for suddenly perturbed 
non-equilibrium systems in m- Although, as compared 
to these studies, we work on systems close to equilibrium 
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Fig. 3. Tpt as a function of integration time t for h = 1 / 4 . The 
physical parameters are Tc = 0.19 and Z\o = 0.29. 



32 64 96 128 160 192 224 256 

t 


Fig. 6. a-s a function of integration time t for h = 

i/i6. The physical parameters are Tc = 0.19 and Ao = 0.083. 



8 16 24 32 40 48 56 64 

t 


Fig. 4. ||at 111 as a function of integration time t for h = i/s. 
The physical parameters are Tc = 0.19 and Aq = 0.16. 



Fig. 7. ipt as a function of integration time t for h = i/ie. The 
physical parameters are Tc = 0.19 and Aq — 0.083. 



8 16 24 32 40 48 56 64 

t 


Fig. 5. 'i/'t as a function of integration time t for h = i/s. The 
physical parameters are Tc = 0.19 and Z\o = 0.16. 


and slightly above the critical temperature, it is inter¬ 
esting to see that our long-term evolutions show oscill- 
lations which resemble the ones predicted for the out-of¬ 
equilibrium case. 


4.4 Results for h = i/ie 

Once more, we depict the time evolution of ||at|| and 
cf. Figs, ini and [T] The conclusions we can draw from these 
two plots are the same as for h = i/s, the only difference 
being the faster oscilllations in line with the bound (1^51) . 
Most importantly, even for this small value of h, we only 
observe diffusion for the linear approximation which, be¬ 
lying its name, does not approximate the BCS equation 
for reasonably long time intervals. Let us summarize our 
results in the concluding Section. 


5 Summary 

We have introduced a reliable integration scheme for the 
time-dependent BCS equation and its linear approxima- 
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tion in spatially homogeneous settings. With the help of 
these algorithms, we could perform numerical long-term 
studies for systems close to equilibrium in order to inves¬ 
tigate the time-evolution of the order parameter at the 
limit close to the critical temperature. The study shows 
very clearly that, opposed to the linear case, the full BCS 
equation does not yield any decay over time in the or¬ 
der parameter Since the conventional time dependent 
Ginzburg-Landau equation is dissipative above the criti¬ 
cal temperature by definition, it cannot give a valid macro¬ 
scopic limit of the full time-dependent BCS/BdG equa¬ 
tions. It can only be seen as the limit of the linearization 
of the full equations but the effects of this linearization 
could clearly be shown not to be negligible in the consid¬ 
ered regime. We thus confirm the analysis provided in [T] . 

In addition, when evolving the system as described by 
the non-linear BGS/BdG equations, we observed oscilla¬ 
tions in the Gooper pair density and in the order param¬ 
eter about a finite value which are similar to oscillations 
which have been observed for out-of-equilibrium systems 
in various works. 

Acknowledgements: We want to thank Ghristian Lu- 
bich for useful remarks and discussions. This work was 
partially funded by the DFG grant GRK 1838. 

Declaration of contribution Gontribution to the orig¬ 
inal manuscript was: G. H. 30%, J. S. 70%. 

Gontribution to the update of the manuscript was: C. H. 
50%, J. S. 50%. 


this into Eq. (TT^ together with the just-determined Aq. 
This yields physically realistic conditions which satisfy the 
energy constraint ([HI). 


B Numerical treatment of the equations 

We want to model a system of infinite spacial extension, 
which, of course, is not possible to achieve on a machine. 
Therefore, we pretend our system to be periodic in space 
but with a large enough period. 


B.l Finite extension and discrete system 

In the Ginzburg-Landau regime, one often takes into ac¬ 
count external potentials that vary on a scale of 0{^/h) 
and, consequently, lead to variations of the system which 
occur over intervals of that very scale. Thus, a valid model 
system should have an extension no smaller than those 
physical variations. But, in order to avoid artificial effects 
due to the periodicity, it is necessary to enlarge this exten¬ 
sions by some multiples of i/zi. For convenience we further¬ 
more include a factor of 27r, wherefore we consider systems 
with period 1 < iV € N. The kernels of the den¬ 
sity operators are now functions on L^([0, i—>■ R). In 

order to simplify the notation, we introduce macroscopic 
variables via Xmac ■= ^jNx. We end up in a 27r-periodic 
setting for which the inner product of two functions / and 
g is just 


A Criticial temperature and initial energy gap 


(27) 

feez 


For translation invariant systems with contact interaction, 
the cricital temperature Tc is well-known to be given im¬ 
plicitly by 



see, e.g. [i®[lD]. The energy gap A between the super¬ 
conducting state and the normal state at temperatures 
beneath the cricital temperature, in turn, can be obtained 
from the relation 


27r 

a 


t,„h 1 ViEpil!) 

~7W^W+W~ 


dp. 


(26) 


In order to calculate the critical temperature and a real¬ 
istic initial value for the gap parameter we thus proceed 
as follows: For a given value of the small crucial parame¬ 
ter h, we first determine the critical temperature Tc and 
set T = Tc — h? ■ For this temperature we then search the 
corresponding gap A following the above definition (1261) 
and set Aq = A a,s its initial state. Finally, as we are 
interested in simulations for temperatures slightly above 
the critical temperature, we put T = Tc + h? and insert 


The self-consistent BGS equations are now given by 

^A(fc) = [i7r,(K),A(fc)], fceZ, (28) 

with the Hamiltonian 

\2[VNh*6ct]k(g-^k‘^)) 

and the Fourier transform Ynh of VNh{-) '■= V (^AO- 
Please note that in the present discrete case the convo¬ 
lution of two summable series and bk has to be under¬ 
stood as 


(a * 


b)k =^ak-jbj. 

jez 


(30) 


B.2 The equations for a delta potential 

For systems on a large torus with a contact interaction ([H]) , 
we can easily see that 


Vnh * ctt 


-a{(j)\a ), 


(31) 
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where (f) is the state given by (j){k) = 1 for all integers 
k. With this, the equations of motion take the convenient 
form 


i'jt{k)=2a {(l)\a)atik) - {4>\a) at{k) 


iat{k) = 2 



at{k) + 2a{4>\a) {2^t{k) 


(32) 
1 ) 

(33) 


for the nonlinear case and 

tdt(fc) = 2 ^t{k) + 2a {(j)\a) {2%{k) - 1) 

(34) 


for the linear case. 

Up to now we are still left with an infinite-dimensional 
system of equations. In order to solve these numerically, 
we have to introduce a suitable finite-dimensional sub¬ 
space. 


corresponding distribution {x) at the discrete points 
Xj = can be computed efficiently with the well-known 
fast Fourier transform (FFT). As we want to resolve phe¬ 
nomena happening on the microscopic scale 0{^/n), we 
choose 


K = M— (38) 

for a large enough integer M. Let us now explain how we 
solve the system of ODE (I5S1) - (I57I) . 


B.4 Solving the system of ordinary differential 
equations 


We first notice that the Hamiltonian Hpt is self-adjoint. 
Thus, the time-evolution of Fj is a unitary transformation 
and, hence, its eigenvalues are preserved. With regard to 
definition m, the eigenvalues can be readily computed as 


B.3 Space discretization 


As the BCS equations are given in their momentum space 
representation, it is most convenient to use the so-called 
Fourier collocation. This means that a 27r-periodic func¬ 
tion f{x) = is approximated by 




tkx 


(35) 


1=-^ 


Ai,2(p) = i ± y (^7t(p) - 0 +\at{pW, (39) 
and we see that the equality 

(^7t(p) “ + l“tb)P = (^7o(p) “ 0 + l«o(p)P 

(40) 


holds. Solving this for yt, we get 


K 


where the coefficients f^{j) are obtained by the discrete 
Fourier transform of the values fj = / (^V^j)) J = ■•■i 

1. Mathematically speaking we work on the subspace spanned 
by the first K eigenfunctions of the Laplacian on [0,27r]. 

As a consequence, the evolution of the system is given by 
the iF-dimensional system of ordinary differential equa¬ 
tions (ODE) 


72- 


lt{p) = 


i - y'Hp) - |dt(p)P for > p, 


* 7 t (k) = 2 a 
K 

- < 

2 - 


(</> 


AI 


K\^K 


r^)df (fc) - (((>"|a 


)df(fc) , (36) 


K 

Y 


- 1 , 




(fc) = 2 


af 


(k) 


K 

- < 

2 - 


,2 

2a{(l>^\a^) ( 27 f (fc) - l) , 
K 


k < 


2 


where we have defined the auxiliary function 


h{p) := 


70 (p) - 2 


|ao(p)7 


(41) 


(42) 


(37) 


The signs in Eq. dUD can be inferred from the initial values 
we use in this work, cf. (US. They are such that %{p) is 
greater than 1/2 for p > p^ and less than or equal to 1/2 
for p < p^. 

Inserting the discrete analogon of Eq. dm into the 
relevant equation of motion dm, we get the nonlinear 
coupled system of equations 


and accordingly for the linear case. From numerical anal¬ 
ysis it is well known that (IHbll yields a very good approx¬ 
imation to 27r-periodic functions with the discretization 
error decreasing rapidly as a function of AT, see, e.g. [H], 
Chapter III.1.3. 

For practical reasons we set K to be an integer power 
of 2 so that for a given af{j), j = —^/ 2 ,..., ^/2 — 1, the 


idf (fc) = 2 af(k) 

± 4a \Jh{k) - |df (/c)|2,^ ^ ^ y 

(43) 

This said, we now present our time integration algorithm. 
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B.5 Time discretization 



1e-05 

Putting it in a formal way, the system we have to integrate 
is given by 


r dy(i) 
J dt 

\y(o) 

with 

= fiy{t)), 

= yo, 

(44) 

U^e-07 

< 

y = (-^/ 2 ) ... 

aK (i^/2 -I))^ e 

(45) 



The right hand side of our initial problem can be written 
as the sum of two terms, 


/(y) = /i(y) + /2(y), 


(46) 


where fi represents the linear part which resembles the 
kinetic part in the Schrddinger-equation and /2 is the non¬ 
linear part. Let r denote a time step and (1>tJ the smooth 
map between y(0) and y(T). Given the special form (1461) 
of the differential equation, one can approximate nu¬ 

merically by 

= (%2,/i o ^t./2 o (yo)- (47) 

This is the well-known Strang splitting. Applying it suc¬ 
cessively yields an approximation to the exact solution at 
times t = nr, n = 1 , 2 ,..., the error of which decreases 
quadratically as a function of the step size r, see, e.g. [^ . 
Chapter 11.5. 

The advantage of the Strang splitting is that <l>r,h can 
be calculated exactly as 

■. (48) 

As for ^r,/ 2 i it has to be approximated due to the nonlin¬ 
earity. For this, we choose a simple Runge-Kutta scheme 
as proposed by |43j whose numerical error is small com¬ 
pared to the error expacted from the splitting Before 
starting the simulations, we still need to fix the mentioned 
discretization parameters r and K. In our case, K itself 
depends on three parameters, cf. Eq. (1^ . As h is the 
semiclassical parameter we want to vary throughout the 
study, we have to choose reasonable values for the remain¬ 
ing quantities M, N and r. We first consider r. 

B.6 Fixing the time discretization parameter 

The step size has to be chosen small enough for both the 
numerical approximation of ^t ,/2 and the Strang splitting 
to give accurate results. For our simulations it turned out 
that reliable results can only be expacted for a step size in¬ 
versely proportional to K. Playing safe we include a small 
factor and set r = o.i/k. As a measure for the time inte¬ 
grator’s accuracy, we consider the discrete analogon of the 
free energy introduced in Eq. m above, which is given by 


32 

t 


Fig. 8. Relative error AFt of the discretized free energy 
against integration time t in semilogarithmic scale for h = i/s. 
The physical parameters are Tc = 0.19 and Ao = 0.16. 


A short calculation yields that this quantity is conserved 
under the exact flow of the corresponding initial value 
problem. Therefore, the reliability of a numerical integra¬ 
tion scheme can be checked by tracking the relative error 
AFt, defined by 


AFT{t) = 


FriF,^) - Ft{F^) 


FT{ri^) 


(50) 


along the numerical evolution. Recurring to a constant 
of motion as a criterion of accurateness is a much applied 
procedure in various computational fields, see, e.g. gms]. 
Following this line of reasoning, we have verified the accu¬ 
racy of our time integrator for every simulation presented 
below. As an example, we show the plot of AFt corre¬ 
sponding to the simulations of Subsection 14.31 in Fig. [5] 


B.7 Fixing the space discretization parameters 

We have seen in the previous Subsection that the time 
step has to be inversely proportional to the dimension of 
the subspace we are approximating our system on. Fur¬ 
thermore, every time step requires a computational ef¬ 
fort which grows linearly with K. Consequently, the com¬ 
plete CPU time for a simulation over a given time inter¬ 
val [0,tend] is quadratic in K. So the dimension of the 
subspace and, thus, the related N and M should be the 
smallest possible. In order to check how small a M we 
can choose without any significant loss of accuracy, we fix 
N = 8 and h = '^/4 and calculate the cricital temperature 
via the discretized version of Eq. (ESI, 



1 

4” 7^ / f^(^mac)|c^ (^mac)| dXmac 

27r Jq 


2-k 

a 


K/2-I 

E 


k=-s :/’2 


tanh 





— fj, 


(51) 


TS{F^) 

(49) 


for different values of M. The result can be seen in Fig. [51 
For different values of N and h we get the same plot. 
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Fig. 9. Critical temperature Tc as a function of the number of 
momenta per unit volume M for N = 8 and h — 1 / 4 . 



Fig. 10. tpt as a function of integration time t for N = 4: and 
h — i/s. The physical parameters are Tc = 0.19 and Aq = 0.16. 


We see that for M = 256 the critical temperature is 
still slightly too small. However, when comparing the evo¬ 
lutions obtained with M = 256 to the according ones 
for M = 512, the relevant figures are indistinguighable 
from each another. For the sake of efficiency, we thus fix 
M = 256 for the rest of this work. 

As for the extension of our interval, N, we have to 
choose it large enough so that the solution cannot reach 
the boundaries during the simulation. As, by construction, 
we work with a periodic setting, a solution reaching one 
end of the interval would enter again at the other end, thus 
leading to unphysical interference. As an example of this 
numerical artifact, we consider the case h = ^/s,, N = 4 
and plot the modulus of ijjt in Fig. (TU] We observe oscil¬ 
lations for larger t which should not show up in reality, 
cf. Subsection 14.31 Whenever we encountered such an ar¬ 
tifact, we successively increased N by factors of 2 until 
the artifact vanished. 
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